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Motivated by earlier studies of artificial perceptions of light called phosphenes, we analyze 
traveling wave solutions in a chain of periodically forced coupled nonlinear oscillators model¬ 
ing this phenomenon. We examine the discrete model problem in its co-traveling frame and 
systematically obtain the corresponding traveling waves in one spatial dimension. Direct nu¬ 
merical simulations as well as linear stability analysis are employed to reveal the parameter 
regions where the traveling waves are stable, and these waves are, in turn, connected to the 
standing waves analyzed in earlier work. We also consider a two-dimensional extension of 
the model and demonstrate the robust evolution and stability of planar fronts and annihi¬ 
lation of radial ones. Finally, we show that solutions that initially feature two symmetric 
fronts with bulged centers evolve in qualitative agreement with experimental observations of 
phosphenes. 

Keywords: coupled nonlinear oscillators, discrete model, traveling waves 


I. INTRODUCTION 

Electrical stimulation of the retina can produce artificial perceptions of luminance changes called 
phosphenes, which may also arise in early stages of retinal or visual disease [1]. The induction of 
phosphenes is being used to help restore vision or develop visual aids for patients with severely 
compromised vision ElEl, and an understanding of how phosphenes arise and behave could con¬ 
tribute to such efforts. In a detailed experimental study [ 1 ], Carpenter explored electrically induced 
phosphenes in human subjects. Each subject’s eyes were immersed in a saline bath to which an 
alternating current was applied. When a dark object was passed through a subject’s visual field in 
the presence of such stimulation, visual perceptions of line phosphenes occurred in its wake. The 
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lines were observed to move and interact but never cross. This work suggested the presence of a 
bistability of activity states in the system, with the moving lines representing boundaries between 
sets of cells in different activity regimes. 

Drover and Ermentrout [5] developed a one-dimensional model providing a simple representation 
of the phosphenes in Carpenter’s experiments and their motion. In their work, a chain of excitable 
neurons was driven by a spatially uniform periodic stimulus at a frequency higher than the cells 
could follow, inducing a 1 : 2 phase locking with the stimulus. In the absence of coupling, neurons 
could fire on even or odd cycles of the stimulus, resulting in an intrinsic bistability for the forced 
system. Sufficiently strong coupling, even if directionally unbiased, elicited unidirectional traveling 
waves in which cells were recruited to switch phase. Large amplitude forcing of an excitable system 
is a difficult problem to tackle analytically. In view of that difficulty, the more recent analysis of [6] 
assumed that each neuron is intrinsically oscillatory and characterized by a state evolving at half 
the frequency of an applied forcing signal. Using multiple time scale expansion and the Fredholm 
alternative, the authors of [6] derived a more analytically tractable effective model for the time 
evolution of the neurons’ phases. This reduced model, which is quite general and particularly 
interesting in its own right, will be the focus of the present work. 

In [6], a detailed numerical existence and stability analysis was done using XPPAUT [7] for a 
finite chain of coupled phase oscillators. Certain interesting bifurcation phenomena were identified 
including saddle-node and pitchfork bifurcations in a two-dimensional parameter space character¬ 
izing the strength and asymmetry of coupling between nearest-neighbor oscillators. Beyond the 
critical points involved in these bifurcations, direct numerical simulations identified traveling waves 
that are strongly reminiscent of the “recruitment waves” obtained in the original chain of forced 
neural oscillators. It is exactly these traveling waves that we systematically obtain and analyze 
in the present manuscript, by a combination of numerical and, whenever possible, semi-analytical 
techniques. 

Specifically, we consider the exact traveling wave problem, which takes the form of an advance- 
delay differential ordinary equation in the co-traveling frame of such waves and which we introduce 
along with the mathematical formulation of the problem in Section 2. In Section 3, we proceed 
to solve this problem numerically, identifying exact (up to numerical error) traveling and standing 
(zero-speed) waves within the full two-dimensional parameter space used in [6]. We explore the 
stability of these waves in two complementary ways. On the one hand, we consider the traveling 
waves as steady states of the associated advance-delay partial differential equation (PDF). On the 
other hand, we examine them via direct numerical simulations of the original system of ordinary 
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differential equations (ODEs) for the coupled oscillators with the initial condition “distilled” on 
the lattice from the obtained traveling wave solution. In [6], this problem was studied in a limited 
range for the parameter fi measuring the asymmetry of the nearest-neighbor coupling function. 
Here, our analysis extends to all values of /i. We show that traveling waves are stable in certain 
parameter regions that are periodic in ji and are located above a certain curve, below which there 
exist stable standing waves. The stability regions alternate along the direction of fi with regions 
where the waves are unstable due to the instability of the background state as well as a frontal 
instability. Simulations of the lattice system initialized by an unstable traveling wave show that 
the frontal instability results in formation of two fronts that propagate in the opposite directions 
with the same speed as the initial wave. 

In the final part of Section 3, we extend the lattice model to a two-dimensional setting and 
show that the planar fronts obtained from the one-dimensional traveling wave are very robust 
even with a local initial distortion, unlike radial fronts, which are eventually annihilated in the 
dynamical evolution. Finally, we consider the evolution of two symmetric fronts with initially 
bulged centers and show that the resulting dynamics is in agreement with Carpenter’s findings, 
based on observations of phosphenes, that lines form loops instead of crossing through each other 
and that a line does not break apart unless it meets another line. This work adds support to the 
idea that a phase model may provide a useful framework for studying phosphenes and also raises 
interesting mathematical issues about the relation between advance-delay PDF and lattice ODE 
system solutions, and we conclude with a brief discussion of these points in Section 4. 

II. DISCRETE MODEL AND THE TRAVELING WAVE EQUATION 

Consider a chain of N periodically forced oscillators governed by the reduced spatially discrete 
model of stimulated retinal cells derived in [6] : 

0-n — — ^-n) + f{^-n) 

— ^[H{0j-i — 9j)-\-H(9j-^i — 9j)]-\- f (9j)^ for j = —n-h 1,..., V — n — 2 (1) 

9N-n-l = kH{9N-n-2 “ + /(6^W-n-l)- 

Here 9j{t) is the slowly evolving phase of each neuron, f{9) is a 7r-periodic forcing or locking 
function, and H{9) is a 27r-periodic function characterizing the coupling of the nearest neighbors 
and multiplied by the coupling constant A; > 0. The periodicities of f{9) and H{9) represent the 
1 : 2 frequency locking present in the system and can be easily generalized to different types of 
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frequency locking. In [6] the number of oscillators was set to be even, N = 2n, but here we also 
allow it to be odd, = 2n + 1. 

In what follows, we consider prototypical examples of the two functions proposed in [6], 

H{9) = sm{9 + /i) — sin(/i), f{9) = — sin(20), 

where the parameter /i measures the asymmetry of the coupling function H{9). We can identify 
0 with 2tt such that 9j G [0, 2tt) and then the firing of the jth neuron corresponds to 9j crossing 
through some distinguished value, taken for some oscillator models, for example, to be 9j = tt. It 
is not hard to see that in the case of /i = 0, when H{9) — sin 9 is odd, there exist equilibrium split 
state (antiphase) solutions given by 

{9-ri, • • •, O-i, 6>o,..., 9N-n-i) = (0,..., 0, TT,..., tt) and (tt, ..., tt, 0,..., 0). 

If = 2n, these solutions have equal numbers of oscillators firing in each cycle. When ji is nonzero, 
such piecewise constant split states no longer exist. However, there are single-front equilibrium 
split states that are close to the above antiphase state but have boundary layers near the front. In 

[5] , such steady state solutions of Q are obtained numerically by varying /i and k. A prototypical 
example of the resulting diagram in the (/x, k) plane is shown in Fig. adapted from Fig. 4 in 

[6] . There we observe that the primary split state solution is stable in a region of small enough 
k (region A) but becomes unstable through a pitchfork bifurcation at A: = A:*(/i) shown by the 
solid curve. This bifurcation is subcritical for fi above a certain threshold, fi > fic- For fi < fic^ 
the bifurcation is supercritical and gives rise to two other, secondary, stable nonsymmetric split 
states that exist in region B. These states, in turn, disappear through a saddle-node bifurcation 
at A: = A:**(/i) (dashed curve), where A:**(/i) > A:*(/i). When exploring the dynamical byproducts 
of these instabilities, the authors of [6] identified a spontaneous emergence of traveling waves in 
region C. This motivates us to seek (numerically) exact traveling wave solutions for this problem. 

To this end, we consider the infinite chain of oscillators 

9j = k[H{9j_i - 9j) + H{9j+i - 9j)] + f{9j), j G Z. (2) 

We start by seeking a solution of this system in the form 9j{t) = ©( 2 ;,t), where z = j — ct is n 
traveling wave coordinate with wave velocity c and r — such that 0^0 and 0 ^ tt when 
z ^ — oc and oo, respectively, or vice versa. Substitution of this ansatz into Q leads to the partial 
differential advance-delay equation (co-traveling frame PDF) of the form 


0r - c0^ =k[H{Q{z -h 1 ,t) - Q{z,r)) -h H{Q{z - 1,t) - 0 ( 2 ;, r))] -h f{Q{z,T)). (3) 
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FIG. 1; Bifurcation diagram for stationary solutions of the system 0 with N = 2n = 10, as obtained in 0. 
The primary split state equilibria are stable in region A and destabilize via a pitchfork bifurcation (solid 
curve, k = /c*(/i)). At /i < /ic the bifurcation is supercritical and gives rise to a pair of secondary split 
state equilibria that are stable in region B and disappear through a saddle-node bifurcation (dashed curve, 
k = n < He)- 

Traveling wave solutions of 


Ojit) = z = j-ct, (4) 

are stationary solutions of (0 and satisfy the advance-delay ordinary differential equation 

-c(j)'{z) = k[H{(p{z + 1) - 4>{z)) + H{4>{z - 1) - ^(z))] -h f{(piz)). (5) 

Since every translate of the traveling solution is also a solution, the additional pinning condition 
0(0) = 7t/2 is imposed on the nonlinear system in order to fix (i.e., pin) the traveling wave. This 
condition uniquely identifies the solution, as well as the corresponding velocity c. Note that the 
traveling wave ansatz @ implies that 0j{t) — 0j+i(t + 1/c) for any integer j, meaning that the 
phase value 0j{t) is periodic (modulo shifts) with a period of 1/c. 
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III. EXISTENCE, STABILITY AND DYNAMICS 
A. Examples of traveling wave solutions 

We start by considering solutions of the traveling wave Eq. i> satisfying 0(0) = 7r/2, 0(2;) ^ tt 
as 2 : ^ oc, (j)[z) ^ 0 as 2 ; ^ —oc. Eq. © is discretized in the interval [—25, 25] with a forward 
difference approximation of the derivative in the left hand side (for comparison purposes, a centered 
difference scheme was also used). The middle grid point is fixed to satisfy the pinning condition 
0(0) = 7r/2, and the speed c is added as a variable. We solve the resulting system using Newton’s 
method for the traveling wave (j){z) and its velocity c, thus identifying numerically exact (at least 
up to a prescribed tolerance of 10“^^ in the discretized system) solutions. 

Examples of traveling wave solutions are shown in Figs. [^-[^ corresponding to different values 
of the parameters k and /i, characterizing the intersite coupling. In each case, stability of the 
traveling wave solution 0 ( 2 ;) is analyzed by computing the spectrum of the Jacobian matrix obtained 
by linearizing Eq. © about 0 ( 2 :) with a forward difference approximation of the spatial derivative 
(see Appendix for further details); the resulting spectra are shown in Figs. We also check 

stability of the obtained solutions by solving both the ODE system 0 and the advance-delay PDF 
© initialized at the traveling wave. Figs. show space-time evolution of the solution of the 

advance-delay PDF Q with the initial condition 0(2:,O) = 0 ( 2 ;). The solution is obtained with 
the discretization grid used to obtain the traveling wave 0 ( 2 ;). As expected, the results show that 
0 ( 2 ;) is a stationary solution of Eq. ©. Finally, in Figs. @ 1 #. we show the space-time evolution 
of the solution of the ODE system 0 with = 2n-h 1 = 51 oscillators initialized by the traveling 
wave solution 0 ( 2 ;) evaluated at the integer values of 2 ;, 0j{0) = 0(j), j = —25,... ,25. The time 
evolution of the PDE Q and the ODE system Q is accomplished using the classical fourth order 
Runge-Kutta method. 

In Fig. [^, the traveling wave solution <^{z) (solid line) is shown at A: = 2.25 and jji — 0.5, 
which yields c = 0.8123. Fig. indicates that the eigenvalues of the linearization Jacobian of the 
PDE (§, evaluated at this solution, have only negative real parts, implying the linear (spectral) 
stability of the traveling wave solution as a stationary solution of Eq. ©• To further examine 
this prediction, the traveling wave solution is subsequently used as the initial data in Eq. ©, 
0 ( 2 :, 0) = 0(^), and the PDE is solved numerically, with the space-time contour of the evolution 
of 0 ( 2 ;, r) shown in Fig. confirming the robust evolution of this steady state. Finally, Fig. §1 
shows the evolution of the solution of the discrete ODE system 0 initialized by the traveling wave 
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solution, 0j{0) = 0(j). As anticipated by the traveling wave ansatz, the ODE solution is found to 
traverse the lattice at the velocity c = 0.8124, which agrees with the traveling wave velocity very 
well [23] . 



FIG. 2: (a) Traveling wave solution (solid curve) (j){z) of (§ dX k = 2.25 and /jl = 0.5, yielding c = 0.8123, 
and the initial condition (circles) 0j(Qi) = (j){j) on the lattice for the ODE simulations shown in (d). (b) 
Eigenvalues of the Jacobian of the linearization of the co-traveling frame PDE ^ around the stationary 
solution (l){z). (c) The space-time evolution of the solution of PDE ^ with the stationary state as the initial 
condition, 0(^,0) = (d) The space-time evolution of the solution of ODE 0 with initial condition 

^j(O) = (j){j) shown by circles in (a). 


Similar results are shown in Figs. and for the cases oi k — 1.5 and fc = 1.1 at /i = 0.5, with 
traveling wave velocities c = 0.5368 and c = 0.2382, respectively, obtained by solving Eq. These 
are in excellent agreement with the velocities c = 0.5367 and c = 0.2377, respectively, obtained 
from the solution of the ODE system 0 initialized by the traveling wave (d panels). Observe 
that as k decreases at fixed /i, the velocity c of the traveling wave decreases, and the eigenvalues 
shown in Figs. and 0 move closer to the imaginary axis as well. Nonetheless, the obtained 










eigenvalues in all three cases remain in the left half-plane ReA < 0, suggesting that the relevant 
wave is stable as a stationary solution of the co-traveling wave PDE Q. We anticipate that this 
spectral stability of the solution for the PDE implies the dynamical stability of the traveling wave 
as the solution of the ODE system as confirmed by our simulations. However, the general 
demonstration of such a connection is, to the best of our knowledge, an intriguing open problem 
in analysis. 



FIG. 3: Same plots as in Fig. [^but now for the stable traveling wave solution with k = 1.5 and fi = 0.5, 
yielding c = 0.5368. 


B. Existence and stability of traveling waves for the chain 

We now explore more systematically the existence and stability of traveling waves in the (/i, k) 
parameter plane. The existence results for sufficiently small /i (0 < /i < 1.5) are summarized 
in Fig. [^, which shows the contour plot of velocity of the traveling wave computed at different 
/i and k using Eq. and complements the steady state analysis of [6] shown in Fig. For 
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FIG. 4: Same plots as in Fig. [^but now for the stable traveling wave solution with k = 1.1 and fi = 0.5, 
yielding c = 0.2382. 


comparison, pitchfork (yellow) and saddle-node (red) bifurcation curves for a 51-oscillator chain 
(and shown by solid and dashed curves, respectively, in Fig. are also included. To obtain the 
pitchfork curve, stationary solutions of Eq. 0 are obtained using Newton’s method, starting with 
a value k = 0.4 for given fi and using continuation in k until the Jacobian of the relevant system 
becomes singular. To obtain the saddle-node curve, a stable stationary solution is identified in 
region B of Fig. This is obtained by integrating Eq. 0 using an unstable solution initial 

condition over a long time horizon until a stable stationary solution is obtained (as an attractor of 
the dynamics). This solution is used in Newton’s method as a good initial guess in order to obtain 
the stable waveform within region B. Numerical continuation is then used with increasing ji until 
the Jacobian becomes singular. This monoparametric process is repated for different values of k 
within region B. The curves obtained are almost indistinguishable from the ones in Fig. While 
we perform our computations in a finite chain, having examined chains of different sizes, we expect 
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FIG. 5: (a) Contour plot of the velocity c of the traveling waves in the (/i, k) plane obtained by solving 
Eq. (§ and isocontours of c = 0.25 (white), c = 0.5 (magenta), c = 0.8 (black) and c = 1.2 (green), shown 
together with the pitchfork (yellow) and saddle-node (red) bifurcation curves obtained for the equilibrium 
states of the 51-oscillator chain, (b) The curve (open circles) where the traveling waves have (almost) zero 
velocity nearly coincides with the saddle-node (green) and pitchfork (blue) bifurcation curves; parameter 
values used in Figures 2|3|4 are also shown (magenta dots). 


our principal conclusions to persist (essentially without modification) for the case of the infinite 
chain. On the flip side, we can identify traveling waves only in region C, while below it their speed 
degenerates to 0, leading to standing waves (split-state equilibria), in agreement with the findings 
of [6] for the finite chain case (and our discussion of standing wave states above). 

We now explore this comparison in more quantitative detail in Fig. [It- Here the open circles 
mark the curve above which numerical simulations of ODE system Q yield stable traveling waves. 
To obtain these points, we fixed /i and solved for the traveling wave solution (j){z) and its 
velocity at (/i, k) starting with k just above the bifurcation curves and then progressively decreasing 
its value. To verify the velocity of the traveling wave at given (/i, fc), we then used 0j{0) = (/)(j) as 
initial conditions in the simulations of 0 and computed the velocity of the propagating front. The 
open circles were obtained by finding the values of k where the traveling wave solution has zero 
velocity up to the numerical error in both methods. The comparison strongly suggests that the 
disappearance of standing wave solutions of [6] gives rise to the traveling wave solutions analyzed 
here. 

The isocontours associated with different velocities, shown in Fig. [5ji, illustrate how the velocity 
c of traveling waves, as obtained from the PDE, depends on the two parameters fi and k. To 











11 


complement these results with a monoparametric visualization, we also show some plots of c{k^ /i) 
at fixed /i in Fig. [^,b and of c as a function of ji at fixed k in Fig. [^,d. These curves show that 
the velocity increases with k and /i. As the velocity decreases, these curves approach the region 
containing the standing waves. 




(a) = 0.5 (b) /i = 1 




FIG. 6: Velocity c as a function of k at fixed fi (panels a and b) and as a function of /i at fixed k (panels c 
and d). 

We now discuss an important aspect of the obtained spectral pictures (e.g. the ones shown 
in Figs. Et) that concerns stability of the background state with ©o = 0 or ©o = tt. 

Our solutions clearly have support on both of these states, so the spectrum of these homogeneous 
states should be mirrored within that of our (standing or traveling) waves. Substituting 0(z,t) = 
©0 + ev{z^T) in Eq. Q, where e is a small parameter, and linearizing about ©o, one obtains the 
following linear advance-delay partial differential equation with constant coefficients for v(z,r): 


Vr — cvz = kcos{/jL){v{z + 1,t) — 2v{z^t) + v{z — 1,t)) — 2v{z^t). 


( 6 ) 
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FIG. 7: (a) The eigenvalues of the Jacobian for the traveling wave solution at /i = 0.5, k = 1.5 discretized 
on the interval [—200, 200] with n = 2001, 4001 and 16001 nodes, (b) The enlarged version of the n = 2001 
and n = 16001 cases from (a), along with the continuous background spectrum curve given by Eq. , near 
the imaginary axis. 


Seeking solutions of Q in the form v = for real we have 

A = —2{2k cos(/i) sin^(p/2) + 1) + icp. 


(7) 


The equilibrium state Q{z^r) = ©o is unstable when ReA > 0, i.e. when 

kcosin) < 

2 Sim 

for some real p. Since /c > 0, this can happen for example at large enough k and | < /i < ^. 
More precisely, the background solution is always unstable for 



k > -|sec(/i)|, 


cos(/i) < 0. 


( 8 ) 


Additionally, since 0 implies that ReA = —2(2A: cos(/i) sin^(p/2) + 1) and ImA = cp^ these eigen¬ 
values determine the following locus of points in the spectral plane: 


ReA = -2 


|2fc cos(/i) sin^ 




(9) 


which should be traceable in the linearization spectra of the traveling waves considered here. 

For 0 < /X < I and ^ < p < 27r, based on equation Q, the background is always stable, 
so any instability must come from the front itself. For example, Fig. 0 shows the eigenvalues of 
the Jacobian matrix associated with the traveling wave solution at /x = 0.5 and k = 1.5. The 
eigenvalues are given for n = 2001, 4001 and 16001 nodes in our discretization of the co-traveling 
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wave problem (© on the interval [—200,200], with the step size h = 400/(n — 1). The eigenvalues 
of the linearized equation about the background solutions ©o = 0 and ©o = tt are given by Eq. 0. 
For these parameter values, we theoretically have —7.27 < ReA < —2. In Fig. one can see a band 
of eigenvalues folding towards the imaginary axis as n increases. Zooming in on the region near the 
imaginary axis in Fig. [^, we see that the eigenvalues are converging to the background spectrum 
as the continuum limit is approached, as predicted. In the Appendix, we explain further how the 
forward difference scheme we use results in the particular structure of the continuous spectrum 
of the problem (inducing the parabolic shape observed in Fig. [^) and how a centered difference 
scheme would affect the corresponding spectral and stability picture. Similar results are found for 
the cases of fc = 2.25 and fc = 1.1 (not shown here). 

In the above discussion, we extended the analysis of the standing waves in [6] to traveling 
wave solutions. However, somewhat in line with the considerations in [6], this study until now 
was limited to 0 < /i < 1.5. We now extend our analysis to the entire upper half of the (/i,/c) 
plane. So far, we have obtained stability curves (pitchfork and saddle-node) separating traveling 
and standing waves. Extending these notions to a broader /i-interval in the (/i, k) plane, we obtain 
the bifurcation curves for 0 < /i < 10 as shown in Fig. Along with extended pitchfork and 
saddle-node curves we show the curves above which the background is unstable, according to 
the inequality Q obtained above for the linearized problem. The curve separating standing and 
traveling waves coincides with the saddle-node bifurcation curves where they exist and follows the 
pitchfork bifurcation curve otherwise. 



FIG. 8: Extended pitchfork (blue) and saddle-node (green) bifurcation curves and the background instability 
curve (red). Magenta circles denote the parameter values where the stability of the traveling waves was 
probed via direct numerical computations. 
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To explore the stability of traveling wave solutions where they exist, we divide the parameter 
space into several regions corresponding to different intervals of /i. These are the ji regions / 
([f)7r]),m {['K,^]),IV ([^,27r]),y (pvr, and VI ([^,37r]). In the regions from 
II to V, we choose points (/i, k) and solve the traveling wave Eq. © for (j){z) and velocity c. Using 
the forward difference or in some cases the centered difference approximation yields convergence to 
the relevant solutions. As before, we then use the resulting traveling wave as initial data to solve 
the ODE system 0 in order to check its stability and visualize the resulting dynamics. To obtain 
the initial condition, we evaluate the traveling wave solution at 51 integer points, and then place 
it on a grid of 81 points by padding it with Os on the left and tts on the right: 


0 ,( 0 ) = 


0, j = -40,..., -26 
j = -25,...,25 
TT, j = 26, ...,40. 


The padding introduces a small perturbation of the initial curve, which, when the solution is 
unstable, initiates the instability. The points (/i, fc) = (1.8,0.75), (2.7,1), (27r — 2.7,1), (27r — 
1.8, .75,), (6,1.6) and (6.5,1.6) (magenta circles in Fig. are chosen as representative points from 
the regions //, ///, IV and V. The points (/i, k) — (2.7,1) and (27r — 2.7,1) lie in the region where 
the background state was shown to be unstable, above the red curve in Fig. Perturbed traveling 
waves for these points are shown in Fig. [^and Fig. [irrespectively. 

When (/i, k) — (2.7,1) G //, the traveling wave is predicted to move to the right with speed 
c = 0.2233. In panel (a) of Fig.jT the initial condition is shown. In panel (b), the solution is shown 
at time t = 4, when the perturbation of the unstable solution has caused the relevant dynamical 
instability to be manifested in the 0 0 part of the solution. In panel (c), the solution is shown at 

time {t = 8.5), when the instability of the background is manifested at both ends of the domain. 
Finally, panel (d) gives the space-time contour plot up to time t — 50, clearly illustrating the 
destabilization of the background and where it is initiated. 


In Fig. 10, when (/i. A:) = (27r — 2.7,1) G ///, the front is predicted to move with speed 
c — —0.2233. The results are similar to the ones discussed above and observed in Fig. Observe 
that the traveling wave solutions (l)i{z) at (/i, k) — (2.7,1), shown in Fig. [^, and 02(^) in Fig. 10a 
at (/i, k) — (27r — 2.7,1) exhibit the symmetry 02(^) — tt — 0i(—z), and their velocities are equal 
in absolute value and opposite in sign. This symmetry is discussed below. Importantly, in both 
cases the intervals where the padding is added, and hence slight perturbations are introduced, are 
exactly where the instability of the background initially manifests itself (in panels (b) and (c)) 
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(c) 


(d) 


FIG. 9: Unstable traveling wave at (/u, k) = (2.7,1) with c = 0.2233: (a) Initial state {t = 0) obtained 

from the solution of Eq. (§. (b) Solution of Eq. 0 at t = 4 where the left side of the solution, with 

Oj « 0, manifests destabilization, (c) Solution at t = 8.5 where the right side, with Oj ^ tt, also shows 

destabilization, (d) Contour plot of space-time evolution until time t = 50. Here the color code represents 

Oj values on the real line, rather than in [0, 27r) mod 27r. 


numerically. 

Fig.[^ shows time snapshots of the solution of the ODE system 0 initialized by the traveling 
wave (solid blue curve) obtained from Eq. Q at (/i,/c)= (1.8,0.75), which lies in the region II 
below the background instability curve. The space-time contour plot of the solution is shown in 
Fig. [TTJd. The traveling wave is initially traveling to the right at the positive predicted velocity 
c = 0.5493. However, a frontal instability of the traveling wave leads to the formation of two 
fronts near t = 7.28 that can be seen at t = 14.56 (dashed-dotted green curve) and t = 21.85 
(dotted magenta curve) in Fig. 11 r, as well as in the contour plot shown in Fig. Interestingly, 
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FIG. 10: Unstable traveling wave at (/r, k) = (27r —2.7,1) with c = —0.2233: (a) Initial state (t = 0) obtained 
from the solution of Eq. (§. (b) Solution of Eq. 0 at t = 4 where the right side of the solution, with Oj « tt 
imanifests destabilization, (c) Solution at t = 8.5 where the left side, with Oj ^ 0, suffers a similar effect, 
(d) Contour plot of space-time evolution until time t = 50. Again, the color code represents Oj values on 
the real line, rather than in [0, 27r) mod 27r. 

the two fronts propagate in the opposite directions with the same speed as the initial unstable 
traveling wave. A similar instability is observed at (/i,/c) = (27r — 1.8,0.75) in region III in 
Fig. . However, in this case the predicted velocity of the traveling wave is found to be negative, 
c = —0.5493. As explained below, this pair of solutions also exhibits the symmetry mentioned 
above for the solutions shown in Fig and [T^. In contrast to those solutions, which featured 
background instability, no instability of the background is observed in these two examples, and only 
the frontal instability arises. At A: = 0.75, the frontal instability exists at least for 1.65 < /i < g, 
where q satisfies 0.75 = 0.51 sec(g)|; that is, until fi reaches the background instability region. No 
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FIG. 11: Unstable traveling wave at (/i,/c) = (1.8,0.75) with c = 0.5493: (a) Snapshots of the solution 
of Eq. 0 initialized at the traveling wave at times t = 0 (solid blue curve), 7.28 (approximate onset of 
instability, dashed red), 14.56 (dash-dotted green) and 21.85 (dotted magenta), (b) The contour plot of 
space-time evolution of the solution (with 6> G M) until time t = 50. The traveling wave becomes unstable 
after an initial transient propagation period and splits into two fronts propagating in the opposite directions 
with the same speed as the initial wave. 


instability is observed for traveling waves at /i < 1.65 and k = 0.75 propagated until time t — 500. 
The results of our stability investigations are summarized in Table I. We emphasize that curves 
in the (/i, A;)-plane separating stable and unstable traveling waves sampled above do not coincide 
with the vertical lines — tv 12 and /i = 37r/2 bounding regions II and III where unstable waves 
were found. In particular, our numerical results indicate that the curve separating stable traveling 
waves at small fi from the unstable ones in region II is located in region II slightly to the right 
of the vertical line fi — tv j 2. By the symmetry described below we anticipate a similar stability 
boundary in region III slightly to the left of the line jj. — 3 tv/2. 

Recall that the background is stable in the regions IV and V. To explore the traveling wave 
stability in these regions, we solve Eq. at points (6,1.6) and (6.5,1.6) in regions IV and U, 
respectively. The fronts and their space-time evolution are shown in Fig. [^for (/i. A:) = (6,1.6) 


and Fig. 14 for (/i, k) = (6.5,1.6). At the point (6,1.6) in region /U, the wave is stable and travels 
to the left, in accordance with the negative velocity c = —0.2919 predicted by the traveling wave 


equation; see Fig. 13 At the point (6.5,1.6) in region U, the wave is also stable and travels to the 
right with the positive velocity c = 0.1894, as shown in Fig. 

To extend these results to the whole upper half of the parameter plane, we observe that the 
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FIG. 12: Unstable traveling wave at (/x, fc) = (27r — 1.8,0.75) with c = —0.5493: (a) Snapshots of the 
solution of Eq. 0 initialized at the traveling wave at times t = 0 (solid blue curve), 7.28 (dashed red), 14.56 
(approximate onset of instability, dash-dotted green) and 21.85 (dotted magenta), (b) The contour plot of 
space-time evolution of the solution (with 6> G M) until time t = 50. The traveling wave becomes unstable 
after an initial propagation transient period and splits into two fronts propagating in opposite directions 
with the same speed as the initial wave. 



FIG. 13: Stable traveling wave at (/r, k) = (6,1.6) with c = —0.2919: (a) The traveling wave (l){z) (solid line) 
and the initial condition 0j{0) = (p^j) (circles) for the simulation of the ODE system 0 . (b) The contour 
plot of the space-time evolution of the solution of ODE system 0 until time t = 100. In accordance with 
the predicted negative velocity, the traveling wave propagates to the left. 

right hand side of the traveling wave Eq. is 27r-periodic in /i, implying that its solutions are 
also 27r-periodic in /i. For example, the traveling wave solutions in region / are isomorphic to the 
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FIG. 14: Stable traveling wave at (/i, k) = (6.5,1.6) with c = 0.1894: (a) The traveling wave (solid line) 
and the initial condition 0j{0) = (j){j) (circles) for the simulation of the ODE system 0. (b) The contour 
plot of the space-time evolution of ODE system 0 until time t = 100. In accordance with the predicted 
positive velocity, the traveling wave now propagates to the right. 


TABLE I: Stability and velocities of traveling waves 


point 

region 

stability 

velocity 

(2.7,1) 

II 

unstable (background) 

0.2233 

(27r-2.7,1) 

III 

unstable (background) 

-0.2233 

(1.8,0.75) 

II 

unstable (frontal) 

0.5493 

(27r - 1.8, 0.75) 

III 

unstable (frontal) 

-0.5493 

(6,1.6) 

IV 

stable 

-0.2919 

(6.5,1.6) 

V 

stable 

0.1894 


ones in region V {II with VI etc.). Note also that the traveling (and standing) wave solutions 
are symmetric about jji = tt. To show this, we define new variables c = —c, /i = 27r — /i, z = — 2 ;, 
^{z) — Ti — (p{z). This yields (p'{z) = ^'{z)^ (j){z + 1) = tt — ^{z — 1) and (j){z — 1) = tt — ^{z + 1). 
Substituting this in Eq. ^ for our choice of II{9) and f{9) we obtain, after simplification, 

—c^'{z) — k{sm{^{z + 1) — ^{z) + fl) — sin(/i) + sin(0(z — 1) — ^{z) + — sin(/i)) — sin(20(^)), 

which in fact coincides with Eq. (§. Hence, ^(z) is a solution to the traveling wave equation at 
jl — 2 tt — ji with velocity c = —c. In other words, for any A: > 0, given a traveling wave (j){z) 
with fl — Ti -\-b (for some 6), there is a traveling wave at /i = tt — 6 obtained through the above 
transformations. This traveling wave at /i = tt — 5 has the same speed and spectrum as the wave 
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with /i = 7r + 6 but propagates in the direction opposite to that of the original traveling wave. This 
implies that the backward moving waves in regions III and IV are obtained directly from the waves 
in regions II and /, respectively, with the spectrum (and hence stability properties) preserved. The 
observed symmetry obviously exists not only at tt but at (2m — l)7r for any integer m. Note that 
in accordance with these results, the extended curve shown in Fig. [^repeats periodically once jjl 
reaches 27r and is symmetric about — it and fi = Stt. Since the waves under study are 27r-periodic, 
this analysis can be extended to the entire upper half of the (/i, k) plane. In particular, it explains 
the symmetry observed above for the solution pairs at (2.7,1), (27r — 2.7,1) and at (1.8,0.75), 
(27r-1.8,0.75). 

In summary, we have obtained an analytical expression for the boundary of a region in parameter 
space where traveling waves suffer a background instability. Numerical simulations show that as 
parameters are varied within from the zero-speed boundary of the traveling wave region, traveling 
waves can be stable or can undergo either a background instability (within the region predicted 
analytically) or a frontal instability (between the stable region and the background instability 
region). Due to a symmetry in the model corresponding to the periodicity of H and /, all results 
repeat periodically in /i. 


C. Two-dimensional phase equation dynamics 

We now turn to some explorations of the two-dimensional generalization of our oscillator prob¬ 
lem. The extension of Eq. Q to two dimensions has the form 

01,1 = A:(i7(0i,2 - 01,i) + 77(02,1 - 0i,i)) + /(0i,i) 

— k{H( 6 i ^2 — 0 z,i) + 77(0^+i,i — 0 ^,i) -h 77(0^_i,i — 02 ,i)) + /( 0 z,i), i = 2 ,..., 2 n — 1 , 

01 ,, = A:(77(0i,,+i - 01 ,,) + 77(0i,,_i - 0i,,) + 77 ( 02 ,, - 0i,,)) + j = 2,..., 2n - 1, 

0,,, = A:(77(0,,,+i - 0,,,) + 77(0,,,_i - 0,,,) + 77(0,+i,, - 0,,,) + 77(0,_i,, - 0,,,)) + /(0,,,), 
ij = 2 , ..., 2 n- 1 , 

( 10 ) 

02, 2 n — ^(77(0^,2n —1 — 02, 2 n) “1“ 77(0^^i,2r2 ~ ^i, 2 n) 77(0^_i,2r2 ~ ^i, 2 n)) T /(02, 222)5 

i = 2 ,..., 2 n — 1 , 

^ 2 nJ — ^(77(02r2,j+l “ ^ 2 nj) + 77(02r2,j-l “ ^ 2 nj) + 77(02r2-l,j “ ^ 2 nj)) + f{S 2 nj)^ 

J ~ 2 ,..., 2 n — 1 , 

0222,222 ~ ^(77(0272,222—1 ^ 2 n^ 2 n) T 77(0272—1,222 0222 , 222 )) “1“ f {^ 2 n, 2 n) • 
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We start by solving Eq. on the interval [—25, 25] at the parameter values A; = 1.3 and /i = 0.5 


to obtain a one-dimensional traveling wave (p{z) with velocity c = 0.4155. We then solve Eq. (10) 
for 0ij{t) using the classical fourth order Runge-Kutta method with the initial data 


= { 


0, i<3 

(/)(z-3), 4<z<54 
TT, i > 55, 

thus studying the two-dimensional evolution of an initially planar front. The planar front propa¬ 
gates in the horizontal direction at the velocity of the one-dimensional wave, suggesting that the 
solution 0ij{t) = (/)(i — ct)^ which corresponds to a stable traveling wave in the one-dimensional 
problem, is also stable in the two-dimensional setting. To illustrate this stability, we distort a 


segment of the planar front in the initial condition, as shown in the first panel {t = 0) in Fig. 15 


The resulting evolution is shown in the remaining panels of Fig. where for better visualization, 
the range of the i-axis is shifted to the left in the last three panels. It can be seen that over time 
the system “heals” the perturbation and gradually restores its quasi-one-dimensional planar front 
character, while the solution eventually settles into traveling with the velocity predicted by the 
one-dimensional results. 

Next we consider the evolution of a radial front. In the first four panels of Fig. [16} we show 
the evolution of an initially circular front for A: = 1.3 and ji = 0.5. The initial condition is set to 
= TT for (i, j) within the circle of radius 30 centered at (40,40), and ^z,j(0) = 0 outside of 
the concentric circle of radius 36. The initial value of 9ij for (z, j) between these circles is obtained 
by linear interpolation. The front shrinks and is eventually annihilated (i.e., disappears). The 
bottom panel of Fig. p!6| displays horizontal slices of the solution at the initial and final time steps 
shown above. Notice that in this case, apparently, the initial radial profile of the front is gradually 
deformed to conform more suitably to the square symmetry of the underlying lattice grid. Hence 
it appears that linear (planar) fronts are fairly robust in this system, while radial ones are clearly 
not as robust and eventually disappear. 

As noted in [5], Carpenter [1] observed that experimentally induced phosphenes move according 
to the following rules in two dimensions: 

1. Lines never cross through one another. Rather, they combine to form loops. 

2. A line never breaks apart unless it meets another line. 

To test whether our two-dimensional model captures these features, we now consider simulations 


with two symmetric fronts that initially bulge either outward (Fig. 17 and Fig. 19) or inward 
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FIG. 15: Snapshots of the evolution of a perturbed planar front introduced within a two-dimensional domain 
with /c = 1.3, /i = 0.5at times t = 10, 200, 400, 800, 1200 and 1400. The evolution shows the healing of the 
front and the decay of the associated perturbation that restore the dynamically robust planar front traveling 
with velocity c = 0.4155. 


(Fig. [T^ and Fig. 20). These initial conditions are shown in the first panel {t = 0) of each figure. 
The simulation results shown in Figs. iTpO , with A: = 1.3, /i = 0.5 in Fig. [T^ and Fig.[T8|and k — 2^ 
jji = 0.15 in Fig. [T^ and Fig. are consistent with Carpenter’s observations listed above. We 
can see that the outwardly bulging fronts eventually touch near the edges of the domain and form 
one loop. Meanwhile, the inwardly bulging fronts eventually touch near their centers and the lines 
break into two parts. These features are very much in line with the expectations of [5] (compare 
with their Fig. 3). However, it should also be mentioned that [5] posit a third and final feature, 
namely that: 

3. Neighboring lines show a tendency to move in a similar manner. 


Our simulations of the coupled oscillator model did not reveal such a tendency. Whether the model 
can be improved to reflect this feature is a question that remains to be considered in future studies. 
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FIG. 16: Snapshots of the evolution of a circular front with k = 1.3, /i = 0.5 at times t = 4, 13, 23 and 
45. The front shrinks (and is eventually annihilated) as time evolves. In the bottom panel, a horizontal 
section(jf = 40) of the front at t = 4 is denoted by plus signs and a horizontal section of the front at t = 45 
(of smaller width) is denoted by circles. 

IV. CONCLUSIONS AND FUTURE CHALLENGES 

In the present work we have revisited a generic nonlinear lattice model derived in [6] and 
associated with the dynamics of a forced network of coupled oscillators. A specific motivation for 
studying this system comes from its relevance to phosphenes, artificial perceptions of light arising 
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FIG. 17: Snapshots of two outwardly perturbed fronts at times t = 0, 20, 40, 60, 65 and 80 when k = 1.3, 
/i = 0.5. The fronts meet near the domain edges, form a loop, and shrink. 

in the visual system in which contours, possibly representing boundaries between sets of neurons 
in different activity states or phases, emerge and propagate. We complemented the important 
initial steady state analysis of [6] (see also m) by exploring the possibility of traveling waves 
in the system of ordinary differential equations. This led us to introduce the co-traveling frame 
(advance-delay) PDE and study existence and stability properties of the traveling waves as both 
special periodic modulo shift solutions of the original systems of ODEs and stationary solutions 
of the co-traveling frame PDE. Our results on traveling waves complement the work of [6] on the 
existence and stability of standing waves. We found a curve above which the standing waves become 
unstable. This curve agrees with and extends the bifurcation curves obtained in [6] by analyzing 
the equilibrium states of the ODE system. We showed that the instability of standing waves above 
the curve leads to the emergence of stable traveling waves in some parameter regimes, while at 
other parameter values traveling waves exhibit either frontal or background instability. An analysis 
of the background steady states provided information about the spectrum of these waves. From an 
applications perspective, a two-dimensional collection of oscillators, corresponding to cells in the 
retina or in a layer of visual cortex, is most relevant, and for this reason we also considered some 
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FIG. 18: Snapshots of two inwardly perturbed fronts at times t = 0, 20, 50, 65, 72 and 90 when k = 1.3, 
/i = 0.5. The fronts meet near their centers, separate, and evolve into two (upper and lower) parts. 

prototypical examples of two-dimensional evolution. In particular, we demonstrated that perturbed 
planar fronts can heal and resume their planar form, while radial fronts shrink. Our simulations 
of the system that initially has two fronts with bulging centers are in qualitative agreement with 
Carpenter’s experimental observations [3] (see also the discussion of IS])- 

The present work leads to numerous interesting questions for the further exploration of this 
and related systems. In particular, obtaining an analytical handle on the spectrum of the front in 
the co-traveling wave PDE and connecting this spectrum to the stability properties of the original 
ODEs would be extremely valuable from a theoretical perspective, not only in the context of the 
present setting but also for wide additional classes of lattice dynamical problems bearing traveling 
waves, such as generalized Frenkel-Kontorova (see [8] and references therein), Fermi-Pasta-Ulam 
(see e.g. iMm and the discussion of Chapter 1 in m) or nonlinear Schrodinger systems (see 
e.g. [ME] and the discussion of Chapter 21 in m)- On the numerical front, while we explored 
a few prototypical cases of two-dimensional evolution, a better understanding of the stationary 
and traveling states in two dimensions clearly merits further investigation. Another interesting 
direction is to formulate and study solutions and spectra of a two-component problem that may 
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FIG. 19: Snapshots of two outwardly perturbed fronts at times t = 0, 50, 100, 130, 135 and 145 when /c = 2, 
/i = 0.15. As in Fig. the fronts meet near the domain edges, form a loop, and shrink. 

support traveling pulses. 

On the applications side, once the basic properties of this model system are more fully under¬ 
stood, it may serve as a form of computational test bed for exploring and making predictions about 
visual phenomena evoked by electrical stimulation. In this context, it may be interesting to explore 
the influence of changes in the amplitude or qualitative form of the forcing function on wave prop¬ 
agation [3] and to investigate waves induced by presenting traveling wave stimuli (e.g., HSIIIH]), 
representing objects passing through the visual field, in addition to periodic forcing. Indeed, these 
adjustments may help resolve the model’s current failure to capture the observation that neigh¬ 
boring fronts tend to propagate in a similar manner. At a more fundamental level, the present 
class of models appears to be a (modified) overdamped variant of the widely studied, so-called sine 
lattices (see e.g. [20] and, for a discussion of some of the relevant applications, [21]), hence it would 
be particularly interesting to explore hybrid variants of these models having as special case limits 
the overdamped and the undamped cases previously explored. Some of these issues are currently 
under investigation and will be reported upon in future publications. 
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FIG. 20: Snapshots of two inwardly perturbed fronts at times t = 0, 50, 150, 170, 180 and 200 when /c = 2, 


fi = 0.15. As in Fig. 18, the fronts meet near their centers, separate, and evolve into two (upper and lower) 
parts. 
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Appendix: some details on numerical methods 

To solve the traveling wave Eq. Q, we used a second order forward difference scheme and 
sometimes a second order centered difference scheme to approximate the term (l)\z). For a grid 
point j, these approximations are of the form (—30j+40j+i — 0 j+ 2 )/( 2 Ax) and (0j+i —0j_i)/(2Ax), 
respectively, where Ax is the grid spacing. In some cases, the forward difference approximation 
used in solving Eq. Q did not converge but the centered difference approximation did. Whichever 
approximation was used to solve Eq. ^ for the traveling wave 0(2;), it was checked that the solution 
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of the ODE system 0 with the initial condition 0i(O) = produced results that were consistent 
with the obtained traveling wave solution. 

In the following we analyze the stability of the background state in the discretized Eq. 0 with 
the forward and centered difference approximations described above. First, the forward difference 
approximation is analyzed. In Fig. we see that the eigenvalues of the Jacobian for Eq. 0 
approximate the eigenvalues of the continuum background as the number of points is increased, 
although the full structure of the forward difference spectral locus is more complex. The eigenvalues 
for the background can also be obtained in the context of the forward difference as follows. In this 
case, Eq. 0 obtained by linearizing Eq. 0 about the background equilibrium state ©o = 0 or 
©0 = TT is replaced by 


Vr-C 


-E(J + 2,t)+4E(J + 1,t) -3E(J,t) 
2Ax 


( 11 ) 


k cos(m) [{ViJ + q,T)- 2V{J, t) + V{J - q, r))] - 2V{J, r), 


where V{J^t) is the approximation of v{z^t) at a grid point zj = JAx for integer J and q is an 
integer such that qAx — 1. Seeking solutions in the form Vj(r) = ^ where p is the wave 

number, and solving for A, we find that 


A = c 


■ cos(2pAx) + 4 cos (p Ax) 

2A^ 


+ 2k cos /i(cos(p) — 1) — 2 + zc 


■ sin(2pAx) + 4 sin(pAx) 


2Ax 


( 12 ) 


The real and imaginary parts of these eigenvalues parametrized by p are shown by red circles in 


Fig. at k = 2.25, 1.5 and 1.1 and p = 0.5. For comparison, the eigenvalues of the Jacobian 
associated with the traveling wave solution are shown by blue pluses (recall also Fig. where 
these eigenvalues are shown for the case k — l.b and p = 0.5 for different numbers of nodes in the 
discretization.) To solve 0, 2001 nodes are used in [—25, 25]. The plots show that the two sets of 
eigenvalues are close to each other. 

Expanding Eq. in Taylor series at small Ax, we obtain 

A = 2{fccos(/i)(cos(p) - 1) - 1} + i{cp + cy (Ax)^ + 0((Ax)^)), 


which yields Eq. 0 in the limit Ax ^ 0. The principal part of the error in the real and imaginary 
parts is —c^(Ax)^ and c^(Ax)^ respectively. As the wave number p increases, this error pushes 
the real part of A to — oc and the imaginary part to +oo and — oc, again in line with the observations 

of Fig.0 

It is interesting to explore the spectral properties of the background of the traveling wave using 
a centered difference approximation instead of the forward difference approximation. Eq. 0 is 
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k=1.1 



FIG. 21: Plot of the eigenvalues of the Jacobian associated with the linearization of Eq. @ about the 
traveling wave solution (blue pluses), solved on [—25, 25] using 2001 nodes and a forward difference scheme. 


and the eigenvalues for the background equilibrium state (red circles) given by Eq. (12). Here fi = 0.5 and 
k = 2.25, 1.5 and 1.1. 


then replaced by 

V - = k cos(/x)(r( J + g, r) - 2V{J, r) + V{J - q, r)) - 2ViJ, r), 

and the eigenvalues A are given by 

A = 2A:cos(/x)(cos(p)-l)+ (13) 

In this case there is no error in the real part of A. 

Generally, this centered difference approximation is numerically unstable for the advection equa¬ 
tion given by © ; see m Here, we find that the nonlinear term does stabilize it for large enough 
values of A:, but the instability in the numerical method is observed for smaller values of k even 
though the solutions of the traveling wave equation are stable, according to the forward difference 
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FIG. 22: Plot of the eigenvalues of the Jacobian associated with the linearization of Eq. (§ about the travel¬ 
ing wave solution (blue pluses), solved on [—25, 25] using 2001 nodes and centered difference approximation. 
Here /i = 0.5 and k = 2.25, 1.5 and 1.1. 


approximation. For example, in Fig. [^ the eigenvalues obtained using the centered difference ap¬ 
proximation of Eq. Q are presented for fc = 1.1, 1.5 and 2.25 with /i = 0.5. The real parts of the 
eigenvalues mostly lie between —5.86 and —2 for k = 1.1, between —7.27 and —2 for k — 1.5 and 
between —9.90 and —2 for k = 2.25, which agree almost exactly with the continuum background 
theory based on Eq. 0 . According to Fig. however, with fixed /i = 0.5, the centered difference 
approximation predicts that traveling waves become unstable somewhere between k = 2.25 and 
k = 1.5, as some eigenvalues emerge with positive real part due to the instabilities associated with 
the centered difference approximation. In the results presented in the manuscript, care has been 
taken to avoid such spurious instabilities induced by the numerical scheme. 
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